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Abstract 

Background: Alternative splicing (AS) of precursor mRNA (pre-mRNA) is an important gene regulation process that 
potentially regulates many physiological processes in plants, including the response to abiotic stresses such as salt stress. 

Results: To analyze global changes in AS under salt stress, we obtained high-coverage (-200 times) RNA sequencing data 
from Arabidopsis thaliana seedlings that were treated with different concentrations of NaCI. We detected that -49% of 
all intron-containing genes were alternatively spliced under salt stress, 10% of which experienced significant differential 
alternative splicing (DAS). Furthermore, AS increased significantly under salt stress compared with under unstressed 
conditions. We demonstrated that most DAS genes were not differentially regulated by salt stress, suggesting that 
AS may represent an independent layer of gene regulation in response to stress. Our analysis of functional categories 
suggested that DAS genes were associated with specific functional pathways, such as the pathways for the responses 
to stresses and RNA splicing. We revealed that serine/arginine-rich (SR) splicing factors were frequently and specifically 
regulated in AS under salt stresses, suggesting a complex loop in AS regulation for stress adaptation. We also showed 
that alternative splicing site selection (SS) occurred most frequently at 4 nucleotides upstream or downstream of the 
dominant sites and that exon skipping tended to link with alternative SS. 

Conclusions: Our study provided a comprehensive view of AS under salt stress and revealed novel insights into the 
potential roles of AS in plant response to salt stress. 
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Background 

High salinity in soil is a major environmental condition 
that adversely affects crop production worldwide. Today, 
roughly 20% of the world's cultivated land and nearly 
half of all irrigated lands are affected by salinity [1]. High 
concentrations of salt in soil lead to ion imbalances and 
hyperosmotic stress in plants. Understanding the mecha- 
nisms of plant responses to salt stress is fundamentally 
important to the study of plant biology and also vital to 
continued development of rational breeding and genetic 
engineering strategies to improve salt tolerance in crop 
plants. Plant's cellular and molecular responses to salt 
stress have been studied intensively [2,3]- Among these 
responses is the dramatic change in the expression of a 
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large number of plant genes, which are regulated at the 
transcriptional as well as the post-transcriptional levels. 

Alternative pre-mRNA splicing (AS) is an important 
mechanism for regulating gene expression and for increas- 
ing transcriptome plasticity and proteome diversity in eu- 
karyotes [4]. AS is involved in many physiological processes 
in plants, including the response to biotic and abiotic 
stresses [5-7]. Although AS of some stress-responsive genes 
has been reported, large-scale or genome-wide studies of 
AS dynamics under salt stress conditions are still relatively 
scarce. Based on the data from Sanger sequencing of full- 
length cDNA libraries from Arabidopsis plants exposed to 
different stresses such as cold, heat, and salt stress, it was 
found that the number of AS events under stress condi- 
tions (particularly low temperature) was significantly 
higher than the number under normal conditions [8]. 
Another study using a whole-genome tiling array in 
Arabidopsis with various stress treatments identified a 
group of AS events that were associated with stress- 
responsive genes and some essential regulatory genes 
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[9]. These and other studies revealed the involvement of 
AS in response to abiotic stress [7]. The methods used 
in these studies, Sanger sequencing and tiling array, 
however suffer from relatively low resolution when 
compared with the recently developed high-throughput 
RNA sequencing (RNA-seq) methods. As a result, some AS 
events, particularly those with lower abundance, may es- 
cape detection. A recent study using high-throughput RNA 
sequencing was conducted with Arabidopsis plants that 
were exposed to various stresses or were at different devel- 
opmental stages and time points in the diurnal cycle [10]. 
That study mainly focused on the complexity of AS rather 
than on a detailed description of the global changes from 
AS under salt stress conditions [10]. 

To investigate the global dynamics of AS under salt 
stress, in this study, we used the Illumina HiSeq platform 
to perform pair-end RNA sequencing with Arabidopsis 
plants that were exposed to different concentrations of salt 
and generated -110 million pair-end reads (101 bp in 
length). In what follows, we first describe the features of AS 
under salt stress based on comparative AS analysis. We 
then report on how the genes with differential AS are well 
associated with specific functional categories, such as the 
responses to stresses and RNA splicing. We also suggest 
that AS could represent a regulatory mechanism inde- 
pendent of the regulation of gene transcriptional activa- 
tion. Finally, we discuss the change in pre-mRNA splicing 
patterns of serine/arginine-rich (SR) splicing factors under 
salt stress. 

Results 

Quality analysis of RNA-seq data 

We used the mRNA-sequencing (RNA-Seq) method to 
acquire whole transcriptomes from both NaCl-treated 
and untreated two-week-old Arabidopsis (ecotype C24) 
seedlings at the single-nucleotide resolution. To detect 
salt-induced AS events precisely, we subjected the seed- 
lings to treatments with different concentrations of NaCl 
(0, 50, 150, or 300 mM). We obtained 110 million se- 
quenced reads (101 bp in length) using the Illumina 
High-Seq sequencing system. On average, nearly 89% of 
these reads could be unambiguously aligned to the 
TAIR10 reference genome sequence (Additional file 1). 
To evaluate the quality of the RNA-seq data, we investi- 
gated the proportion of read alignments in the genome, 
the continuity of reads (3V5 1 bias) along transcriptional 
units (TUs) and sequencing saturation. Firstly, compar- 
ing the mapped reads to the gene annotation revealed 
that about 98% of the reads were from exonic regions, 
whereas only 2% were mapped to intergenic and intronic 
regions (Figure 1A). This was consistent with the quality 
of the Arabidopsis genome assemblies and annotation. 
Secondly, plotting the coverage of reads along each tran- 
script exhibited a uniform distribution with no obvious 



375' bias, reflecting the high quality of the cDNA librar- 
ies (Figure IB). Lastly, we assessed the sequencing satur- 
ation and found that as more reads were obtained, the 
number of newly discovered genes plateaued (Figure 1C), 
suggesting that extensive coverage was achieved. This 
was also supported by plotting the read coverage along 
each chromosome, which showed extensive transcriptional 
activity in the entire genome (Figure ID and Additional 
file 2). To confirm that the comparison of AS was per- 
formed at the same level, we randomly sampled 18 million 
properly paired mapped reads from each RNA-seq library 
for further analysis. 

Identification of AS events 

To identify AS events, we first predicted splice junctions 
using the software TopHat, which was designed to iden- 
tify exon-exon splice junctions. We initially obtained 
433,475 junctions from the four RNA-seq libraries 
(Additional file 3). After filtering the splice junctions 
by two criteria (for details, see Materials and Methods) 
- an overhang size of more than 20 bp and at least two 
reads to support the splice junctions - a junction data 
set of 397,321 confident junctions that we believe to be 
true splice junctions was obtained. Comparison of the 
junctions in this junction data set to the gene annotation 
(TAIR 10) revealed that about 363,383 (91.5%) junctions 
were previously annotated, and that the remaining 33,938 
(8.5%) were novel junctions that had not been annotated in 
the TAIR10 Database (Additional file 3). 

After comparing all the confident junctions to the an- 
notated genes, we identified all the AS events (including 
2275 cassette exons, 6624 alternative 5' splice sites (SSs), 
9654 alternative 3'SSs, 6 mutually exclusive exons, 253 
coordinate cassette exons, 18 alternative first exons and 
10 alternative last exons) under salt stress (Figure 2A). 
We also identified 35,565 intron retention events that had 
at least five intron-reads (i.e., the reads were mapped within 
introns) and more than 80% of the intron region covered 
by intron-reads. Among all these AS events, 45.3% had 
already been annotated in Arabidopsis genes (TAIR10), and 
the remaining 54.7% were identified as novel AS events. 
Based on all identified events, we found that about 49.4% of 
intron-containing genes were alternatively spliced under 
salt stress. 

Intron retention was the most prevalent AS event 
under salt stress, although most intron retentions had 
relatively low read coverage compared to the read cover- 
age of exons (Figure 2B). This is consistent with the 
intron-retention background in Arabidopsis that was re- 
cently reported [11]. Following intron retention, the al- 
ternative 5' and 3' splice sites were relatively prevalent 
compared with the other types of AS events. Sequence 
analysis of alternative 5' splice sites (5'SSs) and alterna- 
tive 3' splice sites (3'SSs) revealed that these activated 
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splice sites were still associated with GU and AG dinucleo- 
tides (Figure 2C). Moreover, we found that the occurrence 
of these alternative 5'SSs and 3'SSs was enriched in the 
downstream and upstream 4 bp region of the dominant 
5'SSs and 3'SSs (Figure 2D), respectively. These features of 
alternative 5'SSs and 3'SSs are consistent with those found 
in the human genome [12]. It is noteworthy that when cor- 
relating exon skipping events to alternative 5'SSs and 3'SSs, 
we found that about -17% of the skipped exons simultan- 
eously had alternative 5'SSs or 3'SSs. This percentage of oc- 
currence was significantly higher than that expected for 
random sampling of all annotated exons (the probability of 
random occurrence is 0.02%, Fisher Exact Test, p < 0.001). 



This result suggests a coordinated occurrence of exon- 
skipping and alternative splice site selection. 

Salt-stress enhances AS 

We next compared the difference in AS between the 
control and NaCl treatments. We found that the number 
of AS events in salt-treated plants was obviously higher 
than that in the control plants (Figure 3), consistent with a 
previous report [8]. We ran a Fishers Exact Test on the 
junction-read-counts/intron-read-counts (only for intron 
retention) and the corresponding exon-read-counts be- 
tween the control and the treatments and identified 2065 
AS events (including 279 alternative 5'SSs, 486 alternative 
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Figure 1 Quality and features of the RNA-seq datasets obtained in the current study. (A) Distribution of the RNA-seq reads along annotated 
Arabidopsis genomic features. Among reads that unambiguously match the Arabidopsis genome, more than 97% reads match annotated exons. (B) 
Relation between the RNA-seq read coverage and the length of the transcriptional unit. The x-axis indicates the relative length of the transcripts, and 
the y-axis shows the median depth of coverage. (C) Saturation curve for gene detection. Randomly sampled reads are plotted against the expressed 
genes. The x-axis shows the number of the mapped reads and the y-axis displays the number of the expressed genes. (D) Transcription profiles are 
plotted across the Arabidopsis genome. The distribution of the RNA-seq read density along the chromosome length is shown. Each vertical blue bar 
represents log2 of the frequency of reads plotted against the chromosome coordinates. A schematic drawing of the chromosome and its features is 
shown below the read density. Approximate boundaries of Arabidopsis centromeres are depicted in violet. 
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Figure 2 Features of alternative splicing (AS) in the Arabidopsis genome. (A) Annotation of AS events based on all the confident junctions. 
The AS events include cassette exons, alternative 5'SSs, alternative 3'SSs, mutually exclusive exons, coordinate cassette exons, alternative first 
exons and alternative last exons. (B) The distribution of read coverage of the retained introns and the corresponding flanking exons. (C) The 
nucleotide sequences around the alternative 5'SSs and 3'SSs are shown by Weblogo. The results indicate that these alternative 5'SSs and 3'SSs 
were still associated with GU and AG dinucleotides. (D) Distribution of activated alternative 5'SSs and 3'SSs around the dominant ones. These 
alternative 5'SSs/3'SSs are enriched in the downstream or upstream 4 bp region of the dominant 5'SSs/3'SSs (denoted as position 0 on the x-axis). 
5'SS, alternative 5' splice site; 3'SS, alternative 3' splice site. 



3'SSs, 102 exon-skipping, and 1198 intron retention events) 
from 1088 genes that were significantly over-represented in 
NaCl-treated plants (Additional files 4, 5, 6, 7 and 8). In 
contrast, we identified only 1320 AS events (including 184 
alternative 5'SSs, 247 alternative 3'SSs, 53 exon-skipping, 
and 836 intron retention events) from 643 genes that were 
absent from these NaCl-treated plants (Additional files 4, 9, 
10, 11 and 12). These data indicated the overall promotion 
of AS by salt stress. 

Changes in splicing patterns associated with stress response 

To investigate the potential influence of salt-stress-induced 
AS on cellular processes, we analyzed functional categories 
and pathways of the genes with differential AS under salt 
stress. We identified 1636 differential alternative splicing 
(DAS) genes in seedlings treated with 50, 150 or 300 mM 
NaCL of which 28.3% were found in the seedlings from at 
least two of these treatments (Figure 4A). An analysis of 
functional categories using the software DAVID [13,14] re- 
vealed that these differentially spliced genes were involved 
in several biological processes, including responses to abi- 
otic stimulus and RNA processing, suggesting that salt 
stress may impact biological processes through changing 
pre-mRNA splicing (Additional file 13). In particular, the 



response-to-abiotic-stimulus functional category was mark- 
edly increased among the DAS genes, and was observed in 
the seedlings in all the salt stress treatments (Figure 4B and 
Additional file 14). The results suggested that AS under salt 
stress was not a random process. Rather, it was associated 
with the stress response. Indeed, further analysis using 
Mapman [15] suggested that genes with aberrant splicing 
in NaCl-treated seedlings were involved in various stress re- 
sponse pathways, including hormone-signaling pathways, 
MAPK-signaling pathways, and transcription regulation of 
stress responses (Figure 4C and Additional file 15). Notably, 
some important genes (such as ERD10, RD22, ATGSTF10, 
ATCPK32, CIPK3 and ERD14) involved in stress responses 
were differentially alternatively spliced in the NaCl-treated 
plants (Figure 4D). Among them, ATCPK32 is an ABA sig- 
naling component that regulates the ABA-responsive gene 
expression via ABF4 [16], and ERD10 encodes a gene in- 
duced by low temperature and dehydrations [17]. Both 
genes showed decreased retention of their first introns 
under salt stress. In contrast, the other three genes (ERD14, 
RD22, and ATGSTF10) involved in abiotic stress responses 
[18-20] showed increased intron retention under salt 
stress. These intron retention events were validated by 
RT-PCR using intron-flanking primers. The amount of 
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Figure 3 The counts of each type of AS events in the control and 50, 150 or 300 mM NaCI treatments. The number of alternative 5'SSs, 
3'SSs and exon-skipping events are more in the NaCI treatments than in the control treatment. The green/blue bars represent forward and reverse 
sequencing reads. 5'SS, alternative 5' splice site; 3'SS, alternative 3' splice site. 



the corresponding PCR products was either increased 
or decreased under salt stress, consistent with the 
RNA-seq data. 

Sequence analysis of these intron-retained transcripts 
suggested that all of these intron retentions could gener- 
ate pre-mature stop codons. Therefore, the decrease or 
increase in intron retention was predicted to increase or 
decrease the abundance of the functional transcripts, re- 
spectively. Since these genes and many other genes with 
significant intron retentions (Additional files 8 and 12) 
have been suggested to play roles in stress responses and 
they are induced by abiotic stresses, an increase in their 
functional transcript levels (e.g., ATCPK32 and ERD10) is 
likely to have positive effects on salt tolerance, whereas a 
decrease in the functional transcript levels of other genes 
(e.g., ERD14, RD22 and ATGSTF10) could have negative 
effects on salt tolerance in plants. These results suggested 
that alteration of AS in stress-responsive genes might im- 
pact a plant's tolerance to salt stress. 

We further compared the functional categories of DAS 
genes with the functional categories of genes without 
DAS. This comparison clearly revealed that different 
functional categories are over-represented in both popu- 
lations (Additional files 13 and 16). Generally, among genes 
that produce alternative transcripts, several Gene Ontology 
(GO) categories related to stress, such as 'response to metal 
ion', 'response to abiotic stimulus', 'response to cadmium 
ion', 'RNA splicing' and 'RNA processing', are over- 



represented. On the other hand, among genes that are not 
alternatively spliced, several functions related to house- 
keeping, such as 'protein transport', 'DNA repair' and 'cell 
wall organization', are over-represented. The presence of 
genes that undergo AS in the stress-related and RNA pro- 
cessing categories is much higher than the presence of 
genes not alternatively spliced in these categories, fur- 
ther supporting the notion that stress-related genes are 
more predisposed to pre-mRNA processing than are 
genes involved in basic cellular functions (Figure 4E). 

AS and gene expression are separately regulated in 
response to salt stress 

From the RNA-seq data, 1,368, 1,901 and 2,729 genes 
were defined as differentially expressed (DE) in the 50, 
150 or 300 mM NaCI treatments, respectively, relative 
to the control (p < 0.01, or fold change >2 and p < 0.05) 
(Additional file 17). The differentially expressed genes 
identified in our study overlapped with those identified 
by other groups based on Microarray analyses of 
salt-stressed Arabidopsis seedlings (data from the 
Genevestigator database), indicating that the salt-stress- 
induced gene regulation found here was comparable to that 
of other studies (Additional file 18). Interestingly, when 
compared with the DAS genes, only 207 DE genes also ex- 
hibited significantly changed intron retention, alternative 
5'SSs, alternative 3'SSs or exon skipping under NaCI treat- 
ments (Figure 5A). Functional categorization of these 207 



Ding et al. BMC Genomics 2014, 15:431 Page 6 of 14 

http://www.biomedcentral.com/1471-2164/15/431 



50 mM NaCI 150 mM NaCI 
66 



I Gene annotated to 
cciresponding pathway 



Biotic/Abiotie Stress 



468 311 
89 

140 



B 



475 
ZOO mM NaCI 



corresponding gene-term association positively reported 
| corresponding gene-term association not reported yet 




response- to abiotic stimulus 
response to organic substance 
respond to endogenous stimulus, 
response lo inorganic subslanc* 
response to hormone stimulus 
response to metal ion 
response to cadmium ion 
RNA processing 
response to osmolic Slress 
rxotein localization 

ntrogen compound oiosynlna ac process 
response io salt siress 
esLablishment of protem locaUiaLion 
protein transport 

response 10 lemperarlure stimulus 
generation of precursor ma tabod as and < 
response to radiation 
response 10 ItCjN stimulus 
cellular response to stress 
intracellular transport 
organic acid OiOSyntheliC process 
carbonyk: acid btosynthatic process 
secondary snetafcokc process 
response to cold 

response to abscisic acid stimulus 
response to oxidase slress 
mBWA molabolic process 
cellular macromolecula localization 
cellular protein localization 

RNA splicing 




response to metal ion 
response to abiotic stimulus 
response to cadmium ion 
response to inorganic substance 
RNft SpJkirig 
mRNA metabolic process 
mRNA processing 
response to osmotic stress 



response to salt stress 




I without DAS 
I with DAS 



logld [It P value) 



Control 
Treatment 



ERD10(AT1G2MW) 




RD2Z(AT5G25610) 




ATG5TF10IAT2G30870) 



Control 
* Treatment 



ATCPK32(AT3GS7530| 



UBQ3 

Figure 4 (See legend on next page.) 



Control 



Treatment 



Control 
Treatment 



eiPK3(AT2G26980) 



1 / 



ERD1J(AT1G76160J 



Ding et al. BMC Genomics 2014, 15:431 
http://www.biomedcentral.com/1471-2164/15/431 



Page 7 of 14 



r <, 

(See figure on previous page.) 

Figure 4 Genes with abnormal splicing in NaCI treatments are closely associated with stress response and transcriptional activation. 

(A) Venn diagram of the overlap of differential alternative splicing (DAS) genes in 50, 150 or 300 mM NaCI treatments. (B) A two-dimensional 
(2-D) view of the relationship between the genes with abnormal splicing and their functional annotations. The functional classification of genes 
was done by the DAVID software. Top 30 of the functional annotations that were ordered by the enrichment scores were selected for the 2-D 
view. (C) A network generated by Mapman indicates that genes with aberrant splicing are involved in various stress response pathways, including 
hormone-signaling pathways, MAPK-signaling pathways and transcription regulation. (D) The representative AS events in six stress-responsive 
genes validated by RT-PCR and visualized by IGV browser. In the RT-PCR validation, the red asterisk (*) to the right side denotes the alternative 
splice form. The red arrow on the top indicates an increase (pointing upward) or decrease (pointing down) in AS events at that salt concentration. 
In the IGV visualization, the exon-intron structure of each gene is given at the bottom of each panel. The grey peaks above the exon-intron 
structure indicate the RNA-seq read density across the gene. The red arrows represent alternative splice sites. The blue arcs in CIPK3 indicate 
splice junction reads that support the junctions. (E) Enrichment of biological processes in DAS genes and genes without DAS. The top 1 0 functional 
categories in DAS genes are shown. In the stress-related and RNA processing categories, the possibility of genes that undergo AS is much higher than 
the possibility of genes that are not alternatively spliced. 



genes revealed that they are enriched in several important 
functional pathways, such as 'response to abiotic stimulus' 
(Figure 5B), indicating that a subset of genes with critical 
salt stress response functions are regulated both transcrip- 
tionally and post-transcriptionally. 

Nonetheless, these co-regulated genes account only 
for a relatively small portion of all DAS or DE genes in 
Arabidopsis. This finding suggested that AS and gene 
activation could be separately regulated in response to 
salt stress. Indeed, analysis of the DAS and DE genes 
confirmed that the over-represented functional categor- 
ies differed largely between the two groups, revealing 
separate regulation of gene expression and AS in re- 
sponse to salt stress (Additional files 19 and 20). For 



A 



example, some functional categories, e.g., 'RNA splicing' 
and 'RNA processing', were over-represented only 
among DAS genes, while other categories, such as 'tran- 
scription' and 'response to hormone stimulus', were 
found among the DE genes (Figure 5C). 

Frequent alteration of SR splicing factors in AS patterns 
under salt stresses 

Whereas genes involved in RNA splicing are mostly not 
regulated by salt stress at the expression level, these 
genes are frequently alternatively spliced under salt 
stress. AS of these splicing-related genes could therefore 
represent an independent means of regulation of genes 
in response to salt stress. Strikingly, we identified 15 
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Figure 5 _ 
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splicing factors with changes in AS under salt stress 
(Additional file 21). Ten of these splicing factors encoded 
SR (serine/arginine-rich) proteins. SR splicing factors play 
key roles in the execution and regulation of pre-mRNA 
splicing in plants. In Arabidopsis, there are a total of 
18 SR proteins [21,22]. Previous studies suggested that 
pre-mRNAs of SR protein genes were frequently alter- 
natively spliced under environmental stress, which is 
thought to alter the splicing of their targets and result 
in adaptive transcriptome changes in response to en- 
vironmental conditions [23,24]. 

We validated six of all these splicing factors by RT-PCR 
and visualized them using the IGV junction browser 
(Figure 6). Among the visualized genes, four SR genes 
(AT-RSP40, AT-RSP41, AT-RS2Z33 and AT-SCL33) ex- 
hibited a decrease in AS events under salt stress, and 
two SR genes (AT-RSP31 and AT-SCL30A) exhibited an 
increase in AS events under salt stress (Figure 6). The 
intron retentions of AT-RS2Z33 and AT-RSP40 were 
detected in the second intron in plants under the 
control conditions, but they were weakly present in the 
plants treated with NaCl. This was also verified by RT- 
PCR using a forward primer in the third exon and a 
reverse primer in the second exon. The intron reten- 
tion in both genes occurred in the 5'UTR region, which 
would lead to abnormal transcripts with long 5'UTRs 
that could interrupt the translation and lead to reduced 
synthesis of the protein. The decreased intron retention in 
these two genes under salt stress should therefore lead to a 
decrease in the level of these corresponding long 5'UTR 
transcripts, which could consequently increase the abun- 
dance of their functional transcripts and proteins. The in- 
tron retention of genes AT-RSP41 and AT-SCL33 that 
respectively occurred in the third and fourth introns were 
clearly detected in the control, while they were barely 
present in samples treated with 300 mM NaCl treatment. 
This was validated by RT-PCR using intron-flanking 
primers. Sequence analysis revealed that both intron reten- 
tion events introduced premature stop codons (PTCs) and 
generated truncated proteins. Therefore, the decreased in- 
tron retention in both genes under salt stress could lead to 
the decrease in abnormal transcripts and the increase in 
functional transcripts. 

Alternative 3'SS and 5'SS were found in the third in- 
tron of AT-SCL30A (AT3G13570) under the 150 mM 
NaCl treatment, while they were weakly present in the 
control. This observation was validated by RT-PCR using 
a forward primer in the third/fourth exon and a reverse 
one covering the splice junction (Figure 6). Further de- 
tailed analysis revealed that both alternative 3'SS and 
5'SS actually introduced a novel exon (not annotated in 
the TAIR10 Arabidopsis genome) that was inserted into 
the region between the third and fourth exon and gener- 
ated a novel isoform. Sequence analysis for this isoform 



suggested that this exon insertion would generate PTCs 
and thus could encode a truncated protein that was 
composed of 120 amino acids. Therefore, this exon in- 
sertion under salt stress could lead to a decrease in the 
functional transcripts. Nonetheless, it is unclear whether 
this novel isoform has any function. 

Finally, we identified an alternative 3'SS in the second 
intron in AT-RSP31, with an increased level under the 
300 mM NaCl treatment. This observation was validated 
by RT-PCR using a forward primer covering the splice 
junction and a reverse one in the second exon (Figure 6). 
This alternative 3'SS (not annotated in the TAIR10 Ara- 
bidopsis genome) extends the third exon into the next 
intron and thus generates a larger exon. Sequence ana- 
lysis for this isoform suggested that this alternative 3'SS 
would introduce PTCs and thus could encode a trun- 
cated protein. It is also unclear whether this truncated 
protein has any function. 

Discussion 

Through comprehensive transcriptome analysis of high- 
throughput RNA-seq data, in this study, we disclosed 
features of genome-wide AS in Arabidopsis under salt 
stress. Our analysis suggests that 49% of the intron- 
containing genes in Arabidopsis genome are alternatively 
spliced under salt-stress conditions. Moreover, we found 
that AS is increased by salt stress and that 10% of the 
intron-containing genes showed significantly differential 
AS under salt-stress conditions. The analysis of func- 
tional categories demonstrated that genes with differen- 
tial AS are associated with responses to stress and RNA 
splicing. Finally, we observed that genes encoding the 
splicing factors, i.e., SR proteins, are subject to frequent 
and specific AS under salt stress. 

An overview of AS in Arabidopsis under salt stress 

Recent studies using massively parallel RNA sequencing 
revealed that a large percentage of genes in Arabidopsis 
undergo AS [10,11], which potentially could significantly 
increase the plasticity of the transcriptome and prote- 
ome diversity. In this study, we conducted systematic 
analysis of the transcriptome of salt-treated Arabidopsis 
plants. Our data revealed that, under salt stress, 49% of 
the intron-containing genes are alternatively spliced. 
This number is higher than that reported by Filichkin 
et al. (42%) [10], but very close to that reported by Li 
et al. (48%) [25], and lower than a recent report that 
61% of multi-exonic genes were alternatively spliced, as 
determined by a normalized cDNA library that facili- 
tated the detection of AS events in low-abundance tran- 
scripts [11]. This marked AS under salt stress could 
provide molecular plasticity for the plants to adapt to 
stress conditions. 
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In this study, we found that intron retention and alter- 
native 5'SSs/3'SSs are much more prominent than exon 
skipping and other types of AS. These observations are 
consistent with the general view of AS in Arabidopsis 



reported previously [10,11]. Importantly, we uncovered 
two novel features of AS in Arabidopsis. First, alternative 
5'SSs/3'SSs tend to occur around the downstream or up- 
stream 4 bp region of the dominant (conical) 5'SS and 
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3'SS (Figure 2D). A similar AS pattern was also reported 
in the human genome [12], suggesting the conservation 
of this AS pattern in eukaryotes. Second, we found a co- 
ordinated occurrence of exon skipping and alternative 
splice site selection. We thus proposed a model where 
exon skipping and alternative splice site selection are 
coupled. We suggest that all the splice sites surrounding 
the dominant ones have the potential to be used as alter- 
native splice sites. These include the splice sites located 
at the next or last exon, which would thus cause exon 
skipping. Previously, exon skipping and alternative splice 
site selection were usually considered as two independ- 
ent AS events. Few links between them were previously 
reported. The discovery of the linkage between these 
two AS events provides a novel perspective of AS and its 
regulation. 

Are stress-induced changes in splicing patterns stress- 
associated acclimation or damage? 

We found that AS events were obviously increased in 
Arabidopsis under salt stress. This finding is consistent 
with some previous studies on AS under environmental 
stresses [5]. For example, cDNA sequencing results indi- 
cated that the number of AS events was significandy higher 
in Arabidopsis plants exposed to different stresses, particu- 
larly low temperature, than in control plants [8]. This in- 
creased AS under stress conditions raises an important 
question on whether the increase is an acclimation re- 
sponse or merely a consequence of splicing errors caused 
by stress damage. We tend to believe that the increase 
comes from splicing errors based on the following reasons. 

First, in another study on the effects of the depletion 
and overexpression of one core component of the spli- 
cing machinery (SAD1, a Sm-like protein 5) on pre- 
mRNA splicing and stress tolerance [26], we found that 
the increase or decrease of AS in many stress-related 
genes can be dynamically controlled by the dosage of 
SAD1; moreover, the increase and decrease in AS are 
closely linked to the sensitivity and tolerance of the 
plants to stress, respectively. Therefore, we considered 
that increased AS could be a result of inaccurate spli- 
cing, which could weaken the function of the corre- 
sponding genes by decreasing the functional transcripts. In 
contrast, decreased AS could be an acclimation to stress re- 
sistance. Secondly, we did observe a stress-induced deregu- 
lation of splicing machinery. In our study, we noticed the 
down-regulation of U6 snRNAs under salt stress in quanti- 
tative RT-PCR assays (Additional file 22). The U6 snRNA is 
a core component of the spliceosome and is required for its 
assembly and catalytic activity during pre-mRNA splicing 
[27,28]. A decrease in the level of this snRNA would likely 
compromise the assembly of the spliceosome and its cata- 
lytic activity [29]. Thirdly, most stress -induced splicing 
variants may not be translated into functional proteins. 



Similarly, some important genes (such as from ERD14, 
RD22 and ATGSTF10) that are involved in abscisic acid 
(ABA) or salt stress responses show increased intron reten- 
tion under salt stress conditions. These transcripts were 
predicted to generate a pre-mature stop codon that would 
lead to non-functional mRNAs or proteins, although we 
currently could not rule out the possibility that some of 
these truncated proteins may still have certain functions in 
plant salt tolerance. Thus, we suggest that stress-induced 
increase in AS could be ascribed to splicing errors or inac- 
curacies caused by stress. 

Nevertheless, if the increase in AS is merely a non- 
specific consequence of stress damage, a random distri- 
bution would be expected among genes. However, our 
data, along with previous reports, demonstrated that the 
genes associated with stress response tend to be alterna- 
tively spliced under stress conditions (Figure 4B). It is 
known that salt stress or other abiotic stresses can activate 
the expression of a large number of plant stress-responsive 
genes that are not expressed or are expressed at lower levels 
under normal non-stressful conditions [30,31]. With the 
simultaneous production of a large amount of these stress- 
inducible pre-mRNAs, cells would need to immediately 
recruit a significant amount of splicing factors and other 
factors for their co-transcriptional or post-transcriptional 
processing. This imposes a huge burden on the splicing 
machinery and, as a result, a significant portion of these 
transcripts fails to be processed adequately when the spli- 
cing machinery is compromised. 

The discussion so far covers only to the global changes 
in AS under salt stress conditions. It should be noted 
that there are indeed specific cases in which AS plays a 
functional role in regulating the response and tolerance 
of plants to stress. Such cases have been described in the 
last few years [review in [5]. This functional role can also 
be seen in the splicing of several SR proteins as dis- 
cussed below. 

Pre-mRNA spicing of SR genes under stress conditions 

The AS pattern of several SR proteins has been shown 
to change obviously under various abiotic stress condi- 
tions, including temperature stress, high salinity and 
high light irradiation [21,32,33]. In this study, we identi- 
fied one- third of the SR genes (six SR genes from four 
SR families) that showed clear changes in AS under salt 
stress. This number is higher than that reported before 
and is probably attributable to the increased sensitivity 
of the sequencing technology used in the current study. 

Interestingly, we clearly identified four SR genes 
(AT-RS2Z33, AT-RSP40, AT-SCL33 and AT-RSP41) that 
showed decreased intron retention under salt stress 
(Figure 6). Sequence analysis revealed that all the 
splice variants with reduced abundance under salt 
stress were aberrant transcripts with premature stop 
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codons that may not produce functional proteins. A 
decrease in these aberrant transcripts and a simultan- 
eous increase in the functional transcripts in these SR 
genes could be an acclimation response to stress that 
may subsequently help to sustain a positive feedback 
loop to increase the splicing efficiency and the produc- 
tion of functional proteins to combat the stress. Con- 
sistently, our recent study demonstrated that the 
mutations of AT-RSP40 or At-RSP41 led to sensitivity 
to salt stress, which implied the positive role of At- 
Rsp40 or At-Rsp41 in salt stress tolerance, probably via 
regulating the pre-mRNA splicing of certain stress tol- 
erance genes [34]. We predict that regulating the ex- 
pression of some of these SR genes or other splicing 
factors may increase plant tolerance to salt stress by 
enhancing the correct splicing of salt tolerance genes. 
Our recent study [26] and a few other studies 
showed that over-expression of certain splicing fac- 
tors indeed could increase plant tolerance to salt 
and other stresses [21,32,33]. 

Conclusions 

Through analyzing global changes in AS under salt 
stress, we firstly identified ~49% of all intron-containing 
genes that were alternatively spliced under salt stress, 
10% of which experienced significant differential alterna- 
tive splicing (DAS). We found that most DAS genes 
were not differentially regulated by salt stress, suggesting 
that AS may represent an independent layer of gene 
regulation in response to stress; DAS genes were associ- 
ated with specific functional pathways, such as the path- 
ways for the responses to stresses and RNA splicing. 
Finally, we revealed that serine/arginine-rich (SR) spli- 
cing factors were frequently and specifically regulated in 
AS under salt stresses, suggesting a complex loop in AS 
regulation for stress adaptation. Therefore, our study 
provided a comprehensive view of AS under salt stress 
and revealed novel insights into the potential roles of AS 
in plant response to salt stress. 

Methods 

Plant materials and growth conditions 

Seeds of the Arabidopsis (ecotype C24) plants were 
surface-sterilized with 50% bleach in 0.01% Triton X-100 
and planted on % Murashige and Skoog (MS) medium 
agar plates supplemented with 3% sucrose. After 4-day 
stratification at 4°C, the plates were placed in a chamber 
(Model CU36-L5, Percival Scientific, Perry, IA, USA) 
under 16 h-white light (~75 umol itT 2 s _1 ) and 8 h-dark 
conditions at 21 ± 1°C for germination and seedling 
growth. Twelve days after being incubated at 21 ± 1°C in 
the chamber, twenty whole seedlings were transferred 
from the agar plate onto filter paper (Catalog No. 05- 
714-4, Fisher Scientific, Pittsburgh, PA, USA) saturated 



with 20 ml of 0, 50, 150, or 300 mM NaCl solution in a 
150 x 15 mm petri dish and incubated in the same 
chamber for 3 h before being harvested and frozen in li- 
quid nitrogen for total RNA extraction [35]. 

RNA extraction, library construction and sequencing 

Using the TRIzol Reagent (Catalog No. 15596-026, Invitro- 
gen), total RNAs were extracted from seedlings without or 
with salt stress treatment. Polyadenylated RNAs were iso- 
lated using the Oligotex mRNA Midi Kit (Catalog No. 
70042, Qiagen). The RNA-seq libraries were constructed 
using the Illumina Whole Transcriptome Analysis Kit fol- 
lowing the standard protocol (Illumina, HiSeq system) and 
sequenced by the Bioscience Core Facility at KAUST on 
the HiSeq platform to generate high-quality pair-end reads. 

Reads alignment and junction prediction 

TopHat [36] was used to align the reads against the 
Arabidopsis genome sequences and annotated gene models 
downloaded from TAIR10 (http://www.arabidopsis.org/) 
with default parameters. Meanwhile, TopHat was also 
used to predict the splice junctions. Based on the gene an- 
notation information, the splice junctions were classified 
into the known and novel splice junctions. In addition, the 
expressed gene or transcripts were identified by the 
Cufflinks software [37]. 

Determination of the criteria for filtering positive junctions 

In the initial prediction, there were a great number of 
novel junctions that had short overhangs (i.e., fewer than 
20 bp) with the corresponding exons, while most of the 
annotated junctions had larger overhangs, with the en- 
richment at -90 bp (Additional file 23 A). Moreover, the 
novel junctions had relatively low coverage compared to 
the annotated junctions (Additional file 23B). In general, 
the junctions with short overhang size and lower cover- 
age were considered as false positives, which are often 
caused by non-specific or error alignment. Therefore, to 
distinguish between true splice junctions and false posi- 
tives, we assessed the criteria based on simulating data 
on a set of randomly-constituted junctions. To do this, 
we firstly generated a set of 80,000 splice junctions in which 
annotated exons from different chromosomes were ran- 
domly selected and spliced together in silico, and 119,618 
annotated junctions from the gene annotation. Since the 
length of our sequencing reads was 101 bp, the splice junc- 
tion sequences were determined to be 180 bp long (90 nt 
on either side of the splice junction) to ensure a 11 bp over- 
hang of the read mapping from one side of the junction 
onto the other. Alignments to the random splice junctions 
were considered to be false positives, because such junc- 
tions are thought to rarely exist when compared to anno- 
tated junctions. The alignments of the raw RNA-seq reads 
to the random junctions revealed that 99.9% of false 
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positive junctions had overhang sizes with fewer than 
20 bp. In sharp contrast, the alignments of the annotated 
junctions indicated that most (98.6%) of annotated junc- 
tions had larger overhang sizes (Additional file 24 A). In 
addition, we estimated that 56.9% of false positive junctions 
had only one read spanning the junction, while the anno- 
tated junctions had higher read coverage (Additional 
file 24B). To minimize the false positive rate, we required 
that the overhang size had to be greater than 20 bp with at 
least two reads spanning the junctions. Using these criteria, 
we filtered out almost all false positive junctions (Additional 
file 24C). 

Annotation of AS events 

JuncBASE [38] was used for annotating all AS events, in- 
cluding cassette exons, alternative 5'SSs, alternative 3'SSs, 
mutually exclusive exons, coordinate cassette exons, alter- 
native first exons, alternative last exons and intron re- 
tention, which takes as input the genome coordinates 
of all annotated exons and all confidently identified 
splice junctions. 

Global comparison of AS 

The global comparison of AS among the control (0), 50, 
150 or 300 mM NaCl treatments was started by equally 
and randomly re-sampling uniquely mapped reads to 
make sure that the comparison was at the same level. 
The comparison refers to the two facets: the absolute 
number of each type of AS event and the number of 
junction reads that was assigned to each type of AS 
event, because both of them can be used to measure the 
global changes of AS. Meanwhile, Fisher's Exact Tests in 
R (http://www.r-project.org/) were used to identify dif- 
ferential representation of each type of AS event, per- 
formed on the number of junction reads that were 
assigned to each type of AS event. 

The identification of differential AS events 

Fisher's Exact Tests were also used to identify the differ- 
ential representation of each AS event. For alternative 
5'SSs and 3'SSs and exon skipping events, Fisher's Exact 
Tests were performed on the comparison of the junction- 
read counts and the corresponding exon-read counts 
between the control and the 50, 150 or 300 mM NaCl 
treatments. The events with p values less than 0.05 were 
identified as significandy differential events. In addition, we 
considered those AS events that were uniquely identified in 
the control or the 50, 150 or 300 mM NaCl treatments sig- 
nificant if there were at least five junction reads to support 
and if the p value of these events was assigned to equal 
zero. Similarly, for intron retention, Fisher's Exact Tests 
were performed on the intron-read counts and the corre- 
sponding exon-read counts between the control and the 50, 
150 or 300 mM NaCl treatments. The events with p value 



less than 0.001 were identified as significant differential 
events. In addition, we considered the intron retention 
events uniquely identified in the control or the 50, 150 or 
300 mM NaCl treatments significant if there were at least 5 
x sequence coverage to support and if more than 80% of 
intron regions were covered by intron-reads. The p value of 
these events was assigned to equal zero. 

RT-PCR validation 

The selected AS and intron retention events were validated 
by RT-PCR using a set of primers (Additional file 25) that 
were designed based on each AS event. Total RNAs from 
the control, 50, 150 or 300 mM NaCl treated seedlings 
were extracted as described above using the TRI solution, 
treated with DNAase I, and reverse-transcribed to cDNA 
(random priming) by using a standard protocol (Super- 
Script II reverse transcriptase, Invitrogen). 

Availability of supporting data 

We deposited the RNA-seq data in the Sequence Read 
Archive (SRA) database with an accession number 
SRP035234. 

Additional files 



Additional file 1: Mapping result of RNA-seq reads. 

Additional file 2: Transcription profiles in the control, 50 or 150 
NaCl treatments were plotted across the Arabidopsis genome. 

Distribution of the RNA-seq read density along the chromosome length 
is shown. Each vertical blue bar represents log2 of the frequency of reads 
plotted against the chromosome coordinates. A schematic drawing of 
the chromosome and its features is shown below the read density. 
Approximate boundaries of Arabidopsis centromeres are depicted in violet. 

Additional file 3: Total junctions from four RNA-seq libraries. 

Additional file 4: The number of DAS events in 50, 150 or 300 mM 
NaCl treatments. Blue bars indicate the number of DAS events that are 
significantly over-represented in NaCl treatment plants. Red bars indicate 
the number of DAS events that are absent in NaCl treatment plants. 

Additional file 5: List of alternative 5' splice sites that were 
over-represented in the 50, 150 or 300 mM NaCl treatments. 

Additional file 6: List of alternative 3' splice sites that were 
over-represented in the 50, 150 or 300 mM NaCl treatments. 

Additional file 7: List of exon-skipping events that were over- 
represented in the 50 mM, 1 50 mM or 300 mM NaCl treatments. 

Additional file 8: List of genes with intron retention in the 50, 150 
or 300 mM NaCl treatments. 

Additional file 9: List of alternative 5' splice sites that were 
over-represented in the control. 

Additional file 10: List of alternative 3' splice sites that were over- 
represented in the control. 

Additional file 11: List of exon-skipping events that were over- 
represented in the control. 

Additional file 12: List of genes with intron retention in the control. 

Additional file 13: Top 50 of the functional categorization 
(biological process) of DAS genes. 

Additional file 14: A two-dimensional (2-D) view of the relationship 
between the genes with abnormal splicing and their functional 
annotations in 50, 150 or 300 mM NaCl treatments. The functional 
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classification of genes was done by the DAVID software. The top 20 
functional annotations that were ordered by the enrichment scores were 
selected for the 2-D view, which indicates that genes with abnormal splicing 
were strikingly enriched in the response-to-abiotic-stress category. 

Additional file 15: A network generated by Mapman indicates that 
genes with aberrant splicing in the 50, 150 or 300 rnM NaCI 
treatments were involved in various stress response pathways, 
including hormone-signaling pathways, MAPK-signaling pathways 
and transcription regulation. 

Additional file 16: Top 50 of the functional categorization 
(biological process) of genes without DAS. 

Additional file 17: List of differentially expressed genes in the 50, 
150 or 300 rnM NaCI treatments relative to the control. 

Additional file 18: The heatmaps generated by mapping the DE 
genes in the 50, 150 or 300 rnM NaCI treatments to the database of 
microarray data with Genevestigator. The heatmaps indicate that DE 
genes (including the up-regulated and down-regulated genes) in NaCI 
treatments are consistent with the gene profiles derived from the 
microarray data. 

Additional file 19: Top 50 of the functional categorization 
(biological process) of DE genes. 

Additional file 20: Top 50 of the functional categorization 
(biological process) of DAS genes. 

Additional file 21: The list of identified SR genes with DAS in 
this study. 

Additional file 22: U6 snRNA expression levels as determined by 
quantitative RT-PCR. The snRNA level in NaCI treatment plants was 
lower than in the control plants. 

Additional file 23: The distinctive features between known and 
novel splice junctions. (A) The density of the overhang size with exons 
for known and novel splice junctions in each sample. The x-axis indicates 
the size of the overhang with exon and the y-axis indicates the density 
of the sizes. A great number of novel junctions has shorter overhangs 
(i.e., fewer than 20 bp) with the corresponding exons, while most of 
the annotated junctions have larger overhang size, with the enrichment 
at -90 bp. (B) The density of junction-read coverage for known and 
novel junctions. The novel junctions have relatively low coverage com- 
pared to the annotated junctions. 

Additional file 24: The features of false positive (random) and 
annotated junctions. (A) The density of the overhang size of false 
positive and annotated junctions. Most of false positive junctions have 
shorter overhang sizes, while the annotated junctions have larger 
overhang sizes. (B) The density of junction-read coverage of false positives 
and annotated junctions. More than half of false positive junctions have only 
one read spanning the junction, while the annotated junctions have 
higher reads coverage. (C) Distinguishing true junctions from false 
positive alignments. To reduce the number of false positive junctions, 
as determined by randomly generated junctions, we reguired that 
the overhang size must be more than 20 bp with at least two reads 
spanning the junctions. Using both criteria, the false positive junctions 
sharply reduced to very low levels (close to zero). In contrast, the 
annotated junctions show no obvious decrease. 

Additional file 25: The primers used for RT-PCR validation. 
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